Biodepth

Dados

Biomassa

BIOMASSA <- read.csv("biodepth/BiodepthSpeciesBiomass.csv", sep = ";", h=T)

Sweden e Sheffield não tem funções de solo medidas & só tem funcoes medidas pra o ano 3 + portugal só sobra 2 amostras quando tira as especies que nao tenho traits = vou filtrar esses

BIOMASSA <- BIOMASSA %>%
  filter(!(site %in% c("Sweden", "Sheffield", "Portugal"))) %>%
  filter(year == 3)

###Biomassa relativa

BIOMASSA <- BIOMASSA %>%
  group_by(year, site, block, plot) %>%
  mutate(percentage_biomass = 100 * biomass / sum(biomass, na.rm = TRUE))

datatable(BIOMASSA, options= list(deferRender = TRUE,  scrollY = "350px", scrollX=T,   scroller = TRUE,  autoWidth=T), rownames=T)
ggplot(BIOMASSA, aes(x = biomass)) +
  geom_histogram(binwidth = 10, fill = "lightblue", color = "black") +
  labs(title = "Biomassa Biodepth", 
       x = "Biomassa", 
       y = "Frequência") +
  theme_minimal()

ggplot(BIOMASSA, aes(x = percentage_biomass)) +
  geom_histogram(binwidth = 10, fill = "lightblue", color = "black") +
  labs(title = "Biomassa relativa Biodepth", 
       x = "Biomassa relativa", 
       y = "Frequência") +
  theme_minimal()

Passando de biomassa para composição de espécies

COMPOSICAO <- BIOMASSA %>% 
  group_by(year,location,site,block,plot,sr,fr,fgc,mix.nest,mix,grass,forb,legume,funct,GRASS,FORB,LEG) %>% 
  pivot_wider(
    names_from = species,
    values_from = percentage_biomass,
    values_fill = list(percentage_biomass=0))

#Germany

COMPOSICAO_GERMANY <- COMPOSICAO %>% 
  filter(site == "Germany")

#matrix de abundancia germany
abundance_germany <- COMPOSICAO_GERMANY %>%
  ungroup() %>%
  mutate(unique_plot = paste(plot, block, year, site, sep = "_")) %>%
  group_by(unique_plot) %>%
  summarise(across(where(is.numeric), sum, na.rm = TRUE)) %>%
  column_to_rownames(var = "unique_plot") %>% 
  dplyr::select(-year,-location, -sr, -fr,-fgc, -mix.nest, -mix, -biomass)
## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `across(where(is.numeric), sum, na.rm = TRUE)`.
## ℹ In group 1: `unique_plot = "B1P001_A_3_Germany"`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
## 
##   # Previously
##   across(a:b, mean, na.rm = TRUE)
## 
##   # Now
##   across(a:b, \(x) mean(x, na.rm = TRUE))
# Cria a coluna 'riqueza' contando o número de espécies em cada plot
abundance_germany$riqueza <- rowSums(abundance_germany > 0)
abundance_germany_filtrado <- abundance_germany[abundance_germany$riqueza > 3, ]
abundance_germany_filtrado$riqueza <- NULL

#Remover espécies que não aparecem em nenhum plot
abundance_germany_filtrado <- abundance_germany_filtrado[, colSums(abundance_germany_filtrado, na.rm = TRUE) > 0]


abundance_matrix_germany <- as.matrix(abundance_germany_filtrado)

#Funcional Germany

# Funcional Germany
SPECIES_germany <- data.frame(species = colnames(abundance_matrix_germany))

SPECIES_germany$species <- tolower(SPECIES_germany$species)

TRAITS$species <- tolower(TRAITS$species) 
TRAITS_germany <- SPECIES_germany %>%
  inner_join(TRAITS, by = "species")

missing_traits_species_germany <- SPECIES_germany %>%
  anti_join(TRAITS, by = "species")

# Diversidade funcional Germany
rownames(TRAITS_germany) <- TRAITS_germany$species
TRAITS_rn_germany <- TRAITS_germany[,-1]

TRAITS_rn_germany$SLA <- log(TRAITS_rn_germany$SLA)
TRAITS_rn_germany$height <- log(TRAITS_rn_germany$height)

TRAITS_rn_germany <- TRAITS_rn_germany[order(rownames(TRAITS_rn_germany)), ]

# Matriz de distância dos traits
euclid_dis_germany <- vegdist(TRAITS_rn_germany, "euclidean")
trait_dat_germany <- as.matrix(euclid_dis_germany)

colnames(abundance_matrix_germany) <- tolower(trimws(colnames(abundance_matrix_germany)))
abundance_matrix_germany <- abundance_matrix_germany[, order(colnames(abundance_matrix_germany))]
fd <- dbFD(trait_dat_germany, abundance_matrix_germany, scale.RaoQ = T,  stand.FRic=T)
## FRic: Dimensionality reduction was required. The last 23 PCoA axes (out of 26 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.9685899
Functional_germany <- data.frame(raoq = fd$RaoQ, FRic = fd$FRic,FEve = fd$FEve, FDis = fd$FDis, FD = fd$nbsp )

datatable(Functional_germany, options= list(deferRender = TRUE,   scroller = TRUE,  scrollX=T, scrollY = "350px", autoWidth=T), rownames=T)
hist(Functional_germany$FRic)

hist(Functional_germany$raoq)

#Greece

COMPOSICAO_GREECE <- COMPOSICAO %>% 
  filter(site == "Greece")

#matrix de abundancia greece
abundance_greece <- COMPOSICAO_GREECE %>%
  ungroup() %>%
  mutate(unique_plot = paste(plot, block, year, site, sep = "_")) %>%
  group_by(unique_plot) %>%
  summarise(across(where(is.numeric), sum, na.rm = TRUE)) %>%
  column_to_rownames(var = "unique_plot") %>% 
  dplyr::select(-year,-location, -sr, -fr,-fgc, -mix.nest, -mix, -biomass)

# Cria a coluna 'riqueza' contando o número de espécies em cada plot
abundance_greece$riqueza <- rowSums(abundance_greece > 0)
abundance_greece_filtrado <- abundance_greece[abundance_greece$riqueza > 3, ]
abundance_greece_filtrado$riqueza <- NULL

#Remover espécies que não aparecem em nenhum plot
abundance_greece_filtrado <- abundance_greece_filtrado[, colSums(abundance_greece_filtrado, na.rm = TRUE) > 0]

abundance_matrix_greece <- as.matrix(abundance_greece_filtrado)

#Funcional Greece

SPECIES_greece <- data.frame(species = colnames(abundance_matrix_greece))

SPECIES_greece$species <- tolower(SPECIES_greece$species)

TRAITS$species <- tolower(TRAITS$species) 
TRAITS_greece <- SPECIES_greece %>%
  inner_join(TRAITS, by = "species")

missing_traits_species_greece <- SPECIES_greece %>%
  anti_join(TRAITS, by = "species")

# Diversidade funcional Greece
rownames(TRAITS_greece) <- TRAITS_greece$species
TRAITS_rn_greece <- TRAITS_greece[,-1]

TRAITS_rn_greece$SLA <- log(TRAITS_rn_greece$SLA)
TRAITS_rn_greece$height <- log(TRAITS_rn_greece$height)

TRAITS_rn_greece <- TRAITS_rn_greece[order(rownames(TRAITS_rn_greece)), ]
euclid_dis_greece <- vegdist(TRAITS_rn_greece, "euclidean")
trait_dat_greece <- as.matrix(euclid_dis_greece)
colnames(abundance_matrix_greece) <- tolower(trimws(colnames(abundance_matrix_greece)))
abundance_matrix_greece <- abundance_matrix_greece[, order(colnames(abundance_matrix_greece))]

fd_greece <- dbFD(trait_dat_greece, abundance_matrix_greece, scale.RaoQ = T,  stand.FRic=T)
## FRic: Dimensionality reduction was required. The last 7 PCoA axes (out of 10 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.9516844
Functional_greece <- data.frame(raoq = fd_greece$RaoQ, FRic = fd_greece$FRic, FEve = fd_greece$FEve, FDis = fd_greece$FDis, FD = fd_greece$nbsp)

datatable(Functional_greece, options= list(deferRender = TRUE, scroller = TRUE, scrollX = T, scrollY = "350px", autoWidth = T), rownames = T)
hist(Functional_greece$FRic)

hist(Functional_greece$raoq)

#Ireland

COMPOSICAO_IRELAND <- COMPOSICAO %>% 
  filter(site == "Ireland")

#matrix de abundancia ireland
abundance_ireland <- COMPOSICAO_IRELAND %>%
  ungroup() %>%
  mutate(unique_plot = paste(plot, block, year, site, sep = "_")) %>%
  group_by(unique_plot) %>%
  summarise(across(where(is.numeric), sum, na.rm = TRUE)) %>%
  column_to_rownames(var = "unique_plot") %>% 
  dplyr::select(-year,-location, -sr, -fr,-fgc, -mix.nest, -mix, -biomass)

# Cria a coluna 'riqueza' contando o número de espécies em cada plot
abundance_ireland$riqueza <- rowSums(abundance_ireland > 0)
abundance_ireland_filtrado <- abundance_ireland[abundance_ireland$riqueza > 3, ]
abundance_ireland_filtrado$riqueza <- NULL

#Remover espécies que não aparecem em nenhum plot
abundance_ireland_filtrado <- abundance_ireland_filtrado[, colSums(abundance_ireland_filtrado, na.rm = TRUE) > 0]

abundance_matrix_ireland <- as.matrix(abundance_ireland_filtrado)

#Funcional Ireland

SPECIES_ireland <- data.frame(species = colnames(abundance_matrix_ireland))

SPECIES_ireland$species <- tolower(SPECIES_ireland$species)

TRAITS$species <- tolower(TRAITS$species) 
TRAITS_ireland <- SPECIES_ireland %>%
  inner_join(TRAITS, by = "species")

missing_traits_species_ireland <- SPECIES_ireland %>%
  anti_join(TRAITS, by = "species")

# Diversidade funcional Ireland
rownames(TRAITS_ireland) <- TRAITS_ireland$species
TRAITS_rn_ireland <- TRAITS_ireland[,-1]

TRAITS_rn_ireland$SLA <- log(TRAITS_rn_ireland$SLA)
TRAITS_rn_ireland$height <- log(TRAITS_rn_ireland$height)

TRAITS_rn_ireland <- TRAITS_rn_ireland[order(rownames(TRAITS_rn_ireland)), ]
euclid_dis_ireland <- vegdist(TRAITS_rn_ireland, "euclidean")
trait_dat_ireland <- as.matrix(euclid_dis_ireland)
colnames(abundance_matrix_ireland) <- tolower(trimws(colnames(abundance_matrix_ireland)))
abundance_matrix_ireland <- abundance_matrix_ireland[, order(colnames(abundance_matrix_ireland))]

fd_ireland <- dbFD(trait_dat_ireland, abundance_matrix_ireland, scale.RaoQ = T,  stand.FRic=T)
## FRic: Dimensionality reduction was required. The last 7 PCoA axes (out of 10 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.9628097
Functional_ireland <- data.frame(raoq = fd_ireland$RaoQ, FRic = fd_ireland$FRic, FEve = fd_ireland$FEve, FDis = fd_ireland$FDis, FD = fd_ireland$nbsp)

datatable(Functional_ireland, options= list(deferRender = TRUE, scroller = TRUE, scrollX = T, scrollY = "350px", autoWidth = T), rownames = T)
hist(Functional_ireland$FRic)

hist(Functional_ireland$raoq)

#Silwood

COMPOSICAO_SILWOOD <- COMPOSICAO %>% 
  filter(site == "Silwood")

#matrix de abundancia silwood
abundance_silwood <- COMPOSICAO_SILWOOD %>%
  ungroup() %>%
  mutate(unique_plot = paste(plot, block, year, site, sep = "_")) %>%
  group_by(unique_plot) %>%
  summarise(across(where(is.numeric), sum, na.rm = TRUE)) %>%
  column_to_rownames(var = "unique_plot") %>% 
  dplyr::select(-year,-location, -sr, -fr,-fgc, -mix.nest, -mix, -biomass)

# Cria a coluna 'riqueza' contando o número de espécies em cada plot
abundance_silwood$riqueza <- rowSums(abundance_silwood > 0)
abundance_silwood_filtrado <- abundance_silwood[abundance_silwood$riqueza > 3, ]
abundance_silwood_filtrado$riqueza <- NULL

#Remover espécies que não aparecem em nenhum plot
abundance_silwood_filtrado <- abundance_silwood_filtrado[, colSums(abundance_silwood_filtrado, na.rm = TRUE) > 0]

abundance_matrix_silwood <- as.matrix(abundance_silwood_filtrado)

#Funcional Silwood

SPECIES_silwood <- data.frame(species = colnames(abundance_matrix_silwood))

SPECIES_silwood$species <- tolower(SPECIES_silwood$species)

TRAITS$species <- tolower(TRAITS$species) 
TRAITS_silwood <- SPECIES_silwood %>%
  inner_join(TRAITS, by = "species")

missing_traits_species_silwood <- SPECIES_silwood %>%
  anti_join(TRAITS, by = "species")

# Diversidade funcional Silwood
rownames(TRAITS_silwood) <- TRAITS_silwood$species
TRAITS_rn_silwood <- TRAITS_silwood[,-1]

TRAITS_rn_silwood$SLA <- log(TRAITS_rn_silwood$SLA)
TRAITS_rn_silwood$height <- log(TRAITS_rn_silwood$height)

TRAITS_rn_silwood <- TRAITS_rn_silwood[order(rownames(TRAITS_rn_silwood)), ]
euclid_dis_silwood <- vegdist(TRAITS_rn_silwood, "euclidean")
trait_dat_silwood <- as.matrix(euclid_dis_silwood)
colnames(abundance_matrix_silwood) <- tolower(trimws(colnames(abundance_matrix_silwood)))
abundance_matrix_silwood <- abundance_matrix_silwood[, order(colnames(abundance_matrix_silwood))]

fd_silwood <- dbFD(trait_dat_silwood, abundance_matrix_silwood, scale.RaoQ = T,  stand.FRic=T)
## FRic: Dimensionality reduction was required. The last 19 PCoA axes (out of 22 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.9660289
Functional_silwood <- data.frame(raoq = fd_silwood$RaoQ, FRic = fd_silwood$FRic, FEve = fd_silwood$FEve, FDis = fd_silwood$FDis, FD = fd_silwood$nbsp)

datatable(Functional_silwood, options= list(deferRender = TRUE, scroller = TRUE, scrollX = T, scrollY = "350px", autoWidth = T), rownames = T)
hist(Functional_silwood$FRic)

hist(Functional_silwood$raoq)

#Switzerland

COMPOSICAO_SWITZERLAND <- COMPOSICAO %>% 
  filter(site == "Switzerland")

#matrix de abundancia switzerland
abundance_switzerland <- COMPOSICAO_SWITZERLAND %>%
  ungroup() %>%
  mutate(unique_plot = paste(plot, block, year, site, sep = "_")) %>%
  group_by(unique_plot) %>%
  summarise(across(where(is.numeric), sum, na.rm = TRUE)) %>%
  column_to_rownames(var = "unique_plot") %>% 
  dplyr::select(-year,-location, -sr, -fr,-fgc, -mix.nest, -mix, -biomass)

# Cria a coluna 'riqueza' contando o número de espécies em cada plot
abundance_switzerland$riqueza <- rowSums(abundance_switzerland > 0)
abundance_switzerland_filtrado <- abundance_switzerland[abundance_switzerland$riqueza > 3, ]
abundance_switzerland_filtrado$riqueza <- NULL

#Remover espécies que não aparecem em nenhum plot
abundance_switzerland_filtrado <- abundance_switzerland_filtrado[, colSums(abundance_switzerland_filtrado, na.rm = TRUE) > 0]

abundance_matrix_switzerland <- as.matrix(abundance_switzerland_filtrado)

#Funcional Switzerland

SPECIES_switzerland <- data.frame(species = colnames(abundance_matrix_switzerland))

SPECIES_switzerland$species <- tolower(SPECIES_switzerland$species)

TRAITS$species <- tolower(TRAITS$species) 
TRAITS_switzerland <- SPECIES_switzerland %>%
  inner_join(TRAITS, by = "species")

missing_traits_species_switzerland <- SPECIES_switzerland %>%
  anti_join(TRAITS, by = "species")

# Diversidade funcional Switzerland
rownames(TRAITS_switzerland) <- TRAITS_switzerland$species
TRAITS_rn_switzerland <- TRAITS_switzerland[,-1]

TRAITS_rn_switzerland$SLA <- log(TRAITS_rn_switzerland$SLA)
TRAITS_rn_switzerland$height <- log(TRAITS_rn_switzerland$height)

TRAITS_rn_switzerland <- TRAITS_rn_switzerland[order(rownames(TRAITS_rn_switzerland)), ]
euclid_dis_switzerland <- vegdist(TRAITS_rn_switzerland, "euclidean")
trait_dat_switzerland <- as.matrix(euclid_dis_switzerland)
colnames(abundance_matrix_switzerland) <- tolower(trimws(colnames(abundance_matrix_switzerland)))
abundance_matrix_switzerland <- abundance_matrix_switzerland[, order(colnames(abundance_matrix_switzerland))]

fd_switzerland <- dbFD(trait_dat_switzerland, abundance_matrix_switzerland, scale.RaoQ = T,  stand.FRic=T)
## FRic: Dimensionality reduction was required. The last 29 PCoA axes (out of 32 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.9723384
Functional_switzerland <- data.frame(raoq = fd_switzerland$RaoQ, FRic = fd_switzerland$FRic, FEve = fd_switzerland$FEve, FDis = fd_switzerland$FDis, FD = fd_switzerland$nbsp)

datatable(Functional_switzerland, options= list(deferRender = TRUE, scroller = TRUE, scrollX = T, scrollY = "350px", autoWidth = T), rownames = T)
hist(Functional_switzerland$FRic)

hist(Functional_switzerland$raoq)

FUNCIONAL juntar

Abundance_juntar

abundance_filtrado <- bind_rows(
  abundance_germany_filtrado,
  abundance_greece_filtrado,
  abundance_ireland_filtrado,
  abundance_silwood_filtrado,
  abundance_switzerland_filtrado
)

MULTIFUNCIONALIDADE

multi <- cbind(dataset_biodep1, getStdAndMeanFunctions(dataset_biodep1, vars))

Dataset final

dataset_final <- multi %>%
  left_join(Functional, by = c("plot", "site"))

datatable(dataset_final, options= list(deferRender = TRUE,    scroller = TRUE, scrollY = "350px", scrollX=T, autoWidth=T), rownames=T)

Análise

cor_matrix
##           raoq       FRic       FEve      FDis         FD
## raoq 1.0000000  0.2061602  0.3401637 0.8031324  0.1210620
## FRic 0.2061602  1.0000000 -0.1454499 0.1812264  0.8138244
## FEve 0.3401637 -0.1454499  1.0000000 0.2740397 -0.2020417
## FDis 0.8031324  0.1812264  0.2740397 1.0000000  0.2830584
## FD   0.1210620  0.8138244 -0.2020417 0.2830584  1.0000000
ggcorrplot(cor_matrix, 
           method = "circle", 
           type = "lower", 
           lab = TRUE, 
           title = "Índices de Diversidade Funcional")

Mixed model

mod1 <- glmmTMB(log(meanFunction) ~ raoq + FRic + (1|site),data = biodepth[-61,])
summary(mod1)
##  Family: gaussian  ( identity )
## Formula:          log(meanFunction) ~ raoq + FRic + (1 | site)
## Data: biodepth[-61, ]
## 
##      AIC      BIC   logLik deviance df.resid 
##     12.2     25.9     -1.1      2.2      108 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.03722  0.1929  
##  Residual             0.05278  0.2297  
## Number of obs: 113, groups:  site, 5
## 
## Dispersion estimate for gaussian family (sigma^2): 0.0528 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.07048    0.09805 -10.918  < 2e-16 ***
## raoq        -0.21096    0.12735  -1.657   0.0976 .  
## FRic         0.55183    0.11579   4.766 1.88e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod_resid <- simulateResiduals(fittedModel = mod1, plot = FALSE)
plot(mod_resid) 

outliers(mod_resid)
## [1] 75
library(MuMIn)
## Warning: pacote 'MuMIn' foi compilado no R versão 4.4.2
r2_mumin <- r.squaredGLMM(mod1)
print(r2_mumin)
##            R2m       R2c
## [1,] 0.1113543 0.4788564
ggplot(biodepth, aes(x = raoq, y = meanFunction)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "RAOQ",
       x = "raoq", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'

ggplot(biodepth, aes(x = raoq, y = meanFunction, color = as.factor(site))) +
  geom_point() +
  geom_smooth(aes(group = as.factor(site)), method = "lm", se = FALSE) +
  labs(title = "RAOQ",
       x = "raoq", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'

ggplot(biodepth, aes(x = FRic, y = meanFunction)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "FRic",
       x = "FRic", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'

ggplot(biodepth, aes(x = FRic, y = meanFunction, color = as.factor(site))) +
  geom_point() +
  geom_smooth(aes(group = as.factor(site)), method = "lm", se = FALSE) +
  labs(title = "FRic",
       x = "FRic", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'

Ceedar Creek

e097

Avaliar o efeito de diferentes níveis de fertilização com nitrogênio em um experimento de competição de árvores.

9 TRTAMENTOS de adição de nitrogênio:
1: 0g/m²/ano nitrogênio (com micronutrientes)
2: 3g/m²/ano nitrogênio (com micronutrientes)
3: 6g/m²/ano nitrogênio (com micronutrientes)
4: 10g/m²/ano nitrogênio (com micronutrientes)
5: 16g/m²/ano nitrogênio (com micronutrientes)
6: 28g/m²/ano nitrogênio (com micronutrientes)
7: 50g/m²/ano nitrogênio (com micronutrientes)
8: 80g/m²/ano nitrogênio (com micronutrientes)
9: 0g/m²/ano (sem nitrogênio sem micronitrientes)

Dados

Biomassa

biomassa <- read.table("e097/e097_Plant aboveground biomass data.txt", h=T, sep="\t")

Filtrando o ano 1982 porque as funções foram coletadas aqui.

biomassa <- biomassa %>%  
  filter(year=="1982")

Biomassa relativa

biomassa <- biomassa %>%
  group_by(Field, plot.number) %>%
  mutate(percentage_biomass = 100 * Species.Biomass..g.m2. / sum(Species.Biomass..g.m2., na.rm = TRUE))

datatable(biomassa, options= list(deferRender = TRUE,  scrollY = "350px", scrollX=T,   scroller = TRUE,  autoWidth=T), rownames=T)
ggplot(biomassa, aes(x = percentage_biomass)) +
  geom_histogram(binwidth = 10, fill = "lightblue", color = "black") +
  labs(title = "Biomassa relativa CEDAR CREEK e97", 
       x = "Biomassa relativa", 
       y = "Frequência") +
  theme_minimal()

###FUNCIONAL

# Cálculo da diversidade funcional
fd <- dbFD(trait_dat, abundance_matrix_clean, scale.RaoQ = T,  stand.FRic=T)
## FRic: Dimensionality reduction was required. The last 60 PCoA axes (out of 64 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.9811237
Functional <- data.frame(raoq = fd$RaoQ, FRic = fd$FRic,FEve = fd$FEve, FDis = fd$FDis, FD = fd$nbsp )

datatable(Functional, options= list(deferRender = TRUE,  scrollY = "350px", scrollX=T,   scroller = TRUE,  autoWidth=T), rownames=T)

MULTIFUNCIONALIDADE

#Dataset junto das funções
funcoes <- potassioFunc %>%
  full_join(fosforoFunc, by = c("Plot.number", "Field.letter")) %>%
  full_join(calcioFunc, by = c("Plot.number", "Field.letter")) %>%
  full_join(magnesioFunc, by = c("Plot.number", "Field.letter")) %>%
  full_join(nitroFunc, by = c("Plot.number", "Field.letter")) %>%
  full_join(phFunc, by = c("Plot.number", "Field.letter"))
e97_multifuncionalidade <- cbind(datset_97, getStdAndMeanFunctions(datset_97, vars))

Dataset final

dataset_final_e97 <- e97_multifuncionalidade %>%
  left_join(Functional, by = c("Plot.number", "treatment", "Field.letter"))

dataset_final_e97$Field.letter <- as.factor(dataset_final_e97$Field.letter)

datatable(dataset_final_e97, options= list(deferRender = TRUE,  scrollY = "350px", scrollX=T,   scroller = TRUE,  autoWidth=T), rownames=T)

Análise

cor_matrix
##            raoq        FRic        FEve       FDis         FD
## raoq  1.0000000  0.51269350 -0.54577474  0.9780512  0.3161771
## FRic  0.5126935  1.00000000 -0.02032126  0.5265112  0.7695444
## FEve -0.5457747 -0.02032126  1.00000000 -0.4516334 -0.0836186
## FDis  0.9780512  0.52651116 -0.45163341  1.0000000  0.3323751
## FD    0.3161771  0.76954441 -0.08361860  0.3323751  1.0000000
ggcorrplot(cor_matrix, 
           method = "circle", 
           type = "lower", 
           lab = TRUE, 
           title = "Índices de Diversidade Funcional")

Mixed model

mod1 <- glmmTMB(log(meanFunction) ~ raoq + FRic + (1|Field.letter),data = e97[-27,])
summary(mod1)
##  Family: gaussian  ( identity )
## Formula:          log(meanFunction) ~ raoq + FRic + (1 | Field.letter)
## Data: e97[-27, ]
## 
##      AIC      BIC   logLik deviance df.resid 
##   -246.8   -233.4    128.4   -256.8      102 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance Std.Dev.
##  Field.letter (Intercept) 0.011419 0.10686 
##  Residual                 0.004854 0.06967 
## Number of obs: 107, groups:  Field.letter, 2
## 
## Dispersion estimate for gaussian family (sigma^2): 0.00485 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.43451    0.07792  -5.576 2.46e-08 ***
## raoq        -0.09278    0.05170  -1.795   0.0727 .  
## FRic        -0.10342    0.08255  -1.253   0.2103    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ATENÇÃO! NÃO TEM HOMOGENEIDADE DA VAR. - nem com log, nem com sqrt
mod_resid <- simulateResiduals(fittedModel = mod1, plot = FALSE)
plot(mod_resid) 

outliers(mod_resid)
## [1] 26
r2_mimine97 <- r.squaredGLMM(mod1)
r2_mimine97
##             R2m       R2c
## [1,] 0.03300648 0.7115699
ggplot(e97, aes(x = raoq, y = meanFunction)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "RAOQ",
       x = "raoq", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'

ggplot(e97, aes(x = raoq, y = meanFunction, color = as.factor(Field.letter))) +
  geom_point() +
  geom_smooth(aes(group = as.factor(Field.letter)), method = "lm", se = FALSE) +
  labs(title = "RAOQ",
       x = "raoq", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'

ggplot(e97, aes(x = FRic, y = meanFunction)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "FRic",
       x = "FRic", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'

ggplot(e97, aes(x = FRic, y = meanFunction, color = as.factor(Field.letter))) +
  geom_point() +
  geom_smooth(aes(group = as.factor(Field.letter)), method = "lm", se = FALSE) +
  labs(title = "FRic",
       x = "Fric", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'

e120

Biodiversity II, effects od plant biodiversity on population and ecossystem process

  • Determinar como a quantidade de espécies de plantas afeta processos ecológicos nos níveis populacional, de comunidade e ecossistêmico. Como a diversidade de plantas influencia a produção de biomassa, a estabilidade do ecossistema, a incidência de doenças, a diversidade de insetos herbívoros e predadores, e outros parâmetros do solo.

168 parcelas de 9m x 9m.

Tratamentos: Parcelas com 1, 2, 4, 8 ou 16 espécies de plantas, com cerca de 30 réplicas para cada nível de diversidade.

COMPOSIÇÃO DAS ESPECIES

datatable(cover, options= list(deferRender = TRUE,  scrollY = "350px", scrollX=T,   scroller = TRUE,  autoWidth=T), rownames=T)
# Calcular a cobertura relativa por quadrante dentro de cada plot
cover <- cover %>%
  group_by(Year, Plot, Quadrat) %>%
  mutate(percentage_cover = 100 * Abscover / sum(Abscover, na.rm = TRUE))

# Somar a cobertura por espécie dentro de cada plot e ano
cover <- cover %>%
  group_by(Year, Plot, Species) %>%
  summarize(Total_Abscover = sum(percentage_cover, na.rm = TRUE)) %>%
  ungroup()
## `summarise()` has grouped output by 'Year', 'Plot'. You can override using the
## `.groups` argument.
# Padronizar para que a soma total seja 100% dentro de cada plot
cover <- cover %>%
  group_by(Year, Plot) %>%
  mutate(cover_final = 100 * Total_Abscover / sum(Total_Abscover, na.rm = TRUE)) %>%
  ungroup()

ggplot(cover, aes(x = Total_Abscover)) +
geom_histogram(bins = 50, fill = "blue", color = "black")

ggplot(cover, aes(x = cover_final)) +
  geom_histogram(bins = 50, fill = "blue", color = "black")

FUNCIONAL

fd_96 <- dbFD(trait_dat_96, abundance_96_matrix, scale.RaoQ = T,  stand.FRic=T)
## FRic: Dimensionality reduction was required. The last 74 PCoA axes (out of 76 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.8967396
Functional_96 <- data.frame(raoq = fd_96$RaoQ, FRic = fd_96$FRic,FEve = fd_96$FEve, FDis = fd_96$FDis, FD = fd_96$nbsp )

hist(Functional_96$FRic)

hist(Functional_96$raoq)

Funcional 2000

SPECIES_00 <- data.frame(species = colnames(abundance_00_matrix))

SPECIES_00$Species <- tolower(SPECIES_00$species)
TRAITS$Species <- tolower(TRAITS$Species) 
TRAITS_00 <- SPECIES_00 %>%
  inner_join(TRAITS, by = "Species") %>% 
  dplyr::select(Species, Leaf.area.per.leaf.dry.mass, Plant.height.vegetative)

rownames(TRAITS_00) <- TRAITS_00$Species
TRAITS_00_rn <- TRAITS_00[,-1]

TRAITS_00_rn$Leaf.area.per.leaf.dry.mass <- log(TRAITS_00_rn$Leaf.area.per.leaf.dry.mass)

TRAITS_00_rn$Plant.height.vegetative <- log(TRAITS_00_rn$Plant.height.vegetative)

TRAITS_00_rn <- TRAITS_00_rn[order(rownames(TRAITS_00_rn)), ]

# Matriz de distância dos traits para o ano 2000
euclid_dis_00 <- vegdist(TRAITS_00_rn, "euclidean")
trait_dat_00 <- as.matrix(euclid_dis_00)

colnames(abundance_00_matrix) <- tolower(trimws(colnames(abundance_00_matrix)))
abundance_00_matrix <- abundance_00_matrix[, order(colnames(abundance_00_matrix))]

colnames(abundance_00_matrix) <- tolower(trimws(colnames(abundance_00_matrix)))
rownames(trait_dat_00) <- tolower(trimws(rownames(trait_dat_00)))

abundance_00_matrix <- abundance_00_matrix[, order(colnames(abundance_00_matrix))]
trait_dat_00 <- trait_dat_00[order(rownames(trait_dat_00)), order(colnames(trait_dat_00))]


# Cálculo da diversidade funcional para o ano 2000
fd_00 <- dbFD(trait_dat_00, abundance_00_matrix, scale.RaoQ = T,  stand.FRic=T)
## FEVe: Could not be calculated for communities with <3 functionally singular species. 
## FDis: Equals 0 in communities with only one functionally singular species. 
## FRic: To respect s > t, FRic could not be calculated for communities with <3 functionally singular species. 
## FRic: Dimensionality reduction was required. The last 71 PCoA axes (out of 74 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.9722786 
## FDiv: Could not be calculated for communities with <3 functionally singular species.
Functional_00 <- data.frame(raoq = fd_00$RaoQ, FRic = fd_00$FRic, FEve = fd_00$FEve, FDis = fd_00$FDis, FD = fd_00$nbsp)

# Visualização dos histogramas para o ano 2000
hist(Functional_00$FRic)

hist(Functional_00$raoq)

# Combina 
Functional_120 <- bind_rows(Functional_96, Functional_00)

MULTIFUNCIONALIDADE

#quantidades de amostragem tudo diferente nos anos 
#carbon_nitroFunc 1164
#biomassaFunc 1862
#soil_carbonFunc 726
#nitrat_amoniFunc 1194
#soil_nitroFunc 726

data_list <- list(carbon_nitroFunc, biomassaFunc, soil_carbonFunc, nitrat_amoniFunc, soil_nitroFunc)
combined_data <- reduce(data_list, full_join, by = c("Year", "Plot"))

#remover todos os dados ausentes
#so sobra 1996 e 2000
funcoes_120 <- na.omit(combined_data)
e120_multifuncionalidade <- cbind(dataset_120, getStdAndMeanFunctions(dataset_120, vars120))

Dataset COMPLETO

dataset_final_120 <- e120_multifuncionalidade %>%
  left_join(Functional_120, by = c("Plot", "Year")) 

datatable(dataset_final_120, options= list(deferRender = TRUE,  scrollY = "350px", scrollX=T,   scroller = TRUE,  autoWidth=T), rownames=T)

Análise

cor_matrix
##            raoq         FRic         FEve       FDis         FD
## raoq 1.00000000  0.349956627  0.044593951 0.87366424 0.07484463
## FRic 0.34995663  1.000000000 -0.006755115 0.34521485 0.49205517
## FEve 0.04459395 -0.006755115  1.000000000 0.04579425 0.04945408
## FDis 0.87366424  0.345214848  0.045794252 1.00000000 0.28825301
## FD   0.07484463  0.492055174  0.049454082 0.28825301 1.00000000
ggcorrplot(cor_matrix, 
           method = "circle", 
           type = "lower", 
           lab = TRUE, 
           title = "Índices de Diversidade Funcional")

Mixed model

mod1 <- glmmTMB(log(meanFunction) ~ raoq + FRic + (1|Year),data = e120)
summary(mod1)
##  Family: gaussian  ( identity )
## Formula:          log(meanFunction) ~ raoq + FRic + (1 | Year)
## Data: e120
## 
##      AIC      BIC   logLik deviance df.resid 
##   -166.4   -147.4     88.2   -176.4      324 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  Year     (Intercept) 3.135e-12 1.771e-06
##  Residual             3.425e-02 1.851e-01
## Number of obs: 329, groups:  Year, 2
## 
## Dispersion estimate for gaussian family (sigma^2): 0.0343 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.94980    0.01921  -49.45   <2e-16 ***
## raoq         0.26408    0.14447    1.83   0.0676 .  
## FRic        -0.08879    0.06768   -1.31   0.1895    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#ATENÇÃO - NAO TEMOS HOMOGENEIDADE DAS VAR. - nem com sqrt nem com log
mod_resid <- simulateResiduals(fittedModel = mod1, plot = FALSE)
plot(mod_resid) 

outliers(mod_resid)
## [1]   3  20  84 115 253
re120 <- r.squaredGLMM(mod1)
re120
##             R2m        R2c
## [1,] 0.01157499 0.01157499
ggplot(e120, aes(x = raoq, y = meanFunction)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "RAOQ",
       x = "raoq", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(e120, aes(x = raoq, y = meanFunction, color = as.factor(Year))) +
  geom_point() +
  geom_smooth(aes(group = as.factor(Year)), method = "lm", se = FALSE) +
  labs(title = "RAOQ",
       x = "raoq", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(e120, aes(x = FRic, y = meanFunction)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "FRic",
       x = "FRic", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(e120, aes(x = FRic, y = meanFunction, color=as.factor(Year))) +
  geom_point() +
  geom_smooth(aes(group = as.factor(Year)), method = "lm", se = FALSE) +
  labs(title = "FRic",
       x = "FRic", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Jena dBEF

Tratamento 1 - Sem história de solo e sem história de plantas: O solo e a camada de plantas foram removidos até 30 cm de profundidade e substituídos por solo de um campo arável adjacente. Tratamento 2 - Com história de solo, mas sem história de plantas: A camada de plantas foi removida, mas o solo original foi mantido. Tratamento 3 - Com história de solo e de plantas: As parcelas existentes do experimento principal, estabelecidas desde 2002, foram mantidas como controle de longo prazo, mantendo tanto a história do solo quanto a das plantas. Esses tratamentos foram projetados para testar como a história do solo e das plantas influenciam a produtividade vegetal e as interações ecológicas ao longo do tempo, considerando tanto a diversidade de plantas quanto os efeitos do histórico do solo sobre a biomassa e outros fatores ecológicos.

80 plots 3 tratamentos 3 anos

COMPOSIÇÃO DAS ESPÉCIES: cobertura

cobertura_2017 <- read.csv("dbef/PLANT_COVER_2017.csv", sep = ";")
cobertura_2018 <- read.csv("dbef/PLANT_COVER_2018.csv", sep = ";")
cobertura_2019 <- read.csv("dbef/PLANT_COVER_2019.csv", sep = ";")

cobertura_total <- bind_rows(cobertura_2017, cobertura_2018, cobertura_2019)
datatable(cobertura_total, options= list(deferRender = TRUE,  scrollY = "350px", scrollX=T,   scroller = TRUE,  autoWidth=T), rownames=T)
## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

FUNCIONAL

#diversidade funcional
fd_dbef <- dbFD(trait_dat_dbef, abundance_dbef_matrix_clean, scale.RaoQ = T,  stand.FRic=T)
## FEVe: Could not be calculated for communities with <3 functionally singular species. 
## FDis: Equals 0 in communities with only one functionally singular species. 
## FRic: To respect s > t, FRic could not be calculated for communities with <3 functionally singular species. 
## FRic: Dimensionality reduction was required. The last 57 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.9033229 
## FDiv: Could not be calculated for communities with <3 functionally singular species.
Functional_dbef <- data.frame(raoq = fd_dbef$RaoQ, FRic = fd_dbef$FRic,FEve = fd_dbef$FEve, FDis = fd_dbef$FDis, FD = fd_dbef$nbsp )

datatable(Functional_dbef, options= list(deferRender = TRUE,  scrollY = "350px", scrollX=T,   scroller = TRUE,  autoWidth=T), rownames=T)

MULTIFUNCIONALIDADE

# Combinar os Dados com a Base
dados_combinados <- base %>%
  left_join(cobertura_media, by = c("plot", "treatment", "year")) %>%
  left_join(biomass_summary, by = c("plot", "treatment", "year")) %>%
  left_join(atv_micro, by = c("plot", "treatment", "year")) %>%
  left_join(herbivoria_total, by = c("plot", "treatment", "year")) %>%
  left_join(predacao, by = c("plot", "treatment", "year"))
dbef_multifunc <- cbind(dados_combinados, getStdAndMeanFunctions(dados_combinados, func_vars))

DATASET COMPLETAO

dataset_final_JENA <- dbef_multifunc %>%
  left_join(Functional_dbef, by = c("plot", "treatment", "year"))

datatable(dataset_final_JENA, options= list(deferRender = TRUE,  scrollY = "350px", scrollX=T,   scroller = TRUE,  autoWidth=T), rownames=T)

Análise

cor_matrix
##           raoq      FRic      FEve      FDis        FD
## raoq 1.0000000 0.5345951 0.2837215 0.9469518 0.3352497
## FRic 0.5345951 1.0000000 0.1763367 0.4895485 0.7310809
## FEve 0.2837215 0.1763367 1.0000000 0.3207152 0.2401041
## FDis 0.9469518 0.4895485 0.3207152 1.0000000 0.3560987
## FD   0.3352497 0.7310809 0.2401041 0.3560987 1.0000000
ggcorrplot(cor_matrix, 
           method = "circle", 
           type = "lower", 
           lab = TRUE, 
           title = "Índices de Diversidade Funcional")

Mixed model

mod1 <- glmmTMB(meanFunction ~ raoq + FRic + (1|year),data = jena[-364,])
summary(mod1)
##  Family: gaussian  ( identity )
## Formula:          meanFunction ~ raoq + FRic + (1 | year)
## Data: jena[-364, ]
## 
##      AIC      BIC   logLik deviance df.resid 
##  -1289.1  -1269.6    649.6  -1299.1      363 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev.
##  year     (Intercept) 0.0007979 0.02825 
##  Residual             0.0016592 0.04073 
## Number of obs: 368, groups:  year, 3
## 
## Dispersion estimate for gaussian family (sigma^2): 0.00166 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.21583    0.01687  12.794  < 2e-16 ***
## raoq        -0.08669    0.03044  -2.848   0.0044 ** 
## FRic         0.04534    0.01038   4.367 1.26e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rjena <- r.squaredGLMM(mod1)
rjena
##            R2m       R2c
## [1,] 0.0343871 0.3479555
mod_resid <- simulateResiduals(fittedModel = mod1, plot = FALSE)
plot(mod_resid)

outliers(mod_resid)
## integer(0)
ggplot(jena, aes(x = raoq, y = meanFunction)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "RAOQ",
       x = "raoq", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 164 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 164 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(jena, aes(x = raoq, y = meanFunction, color = as.factor(year))) +
  geom_point() +
  geom_smooth(aes(group = as.factor(year)), method = "lm", se = FALSE) +
  labs(title = "RAOQ",
       x = "raoq", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 164 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 164 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(jena, aes(x = raoq, y = meanFunction, color = as.factor(year))) +
  geom_point() +
  geom_smooth(aes(group = as.factor(year)), method = "lm", se = FALSE) +
  labs(title = "RAOQ",
       x = "raoq", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 164 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 164 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(jena, aes(x = FRic, y = meanFunction, color = as.factor(year))) +
  geom_point() +
  geom_smooth(aes(group = as.factor(year)), method = "lm", se = FALSE) +
  labs(title = "FRic",
       x = "FRic", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 351 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 351 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(jena, aes(x = FRic, y = meanFunction)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "FRic",
       x = "FRic", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 351 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 351 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(jena, aes(x = FRic, y = meanFunction, color=as.factor(year))) +
  geom_point() +
  geom_smooth(aes(group = as.factor(year)), method = "lm", se = FALSE) +
  labs(title = "FRic",
       x = "FRic", y = "meanFunction") +
  theme_get()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 351 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 351 rows containing missing values or values outside the scale range
## (`geom_point()`).